Unobserved riverflows in Québec, Canada

The functionalities of ErrorsInVariablesExtremes.jl are illustrated by producing an Extreme Value Analysis of unobserved riverflows of a section of the Chaudière River in Québec, Canada, from 1960 to 2020.

Load the data

Loading the pseudo-observations into a vector of Pseudodata:

filename = "../../../test/data/A2020_Analyse_Historique_QMA_SLSO00003.nc"
pdata = load_discharge_distribution(filename)

6-element Vector{Pseudodata}:
 Pseudodata:
  name: LN24HA
  year: Vector{Int64}[60]
  value:Vector{Distributions.LogNormal{Float64}}[60]

 Pseudodata:
  name: MG24HA
  year: Vector{Int64}[60]
  value:Vector{Distributions.LogNormal{Float64}}[60]

 Pseudodata:
  name: MG24HI
  year: Vector{Int64}[60]
  value:Vector{Distributions.LogNormal{Float64}}[60]

 Pseudodata:
  name: MG24HK
  year: Vector{Int64}[60]
  value:Vector{Distributions.LogNormal{Float64}}[60]

 Pseudodata:
  name: MG24HQ
  year: Vector{Int64}[60]
  value:Vector{Distributions.LogNormal{Float64}}[60]

 Pseudodata:
  name: MG24HS
  year: Vector{Int64}[60]
  value:Vector{Distributions.LogNormal{Float64}}[60]

The vector of years can be extracted:

years = pdata[1].year;

60-element Vector{Int64}:
 1961
 1962
 1963
 1964
 1965
 1966
 1967
 1968
 1969
 1970
    ⋮
 2012
 2013
 2014
 2015
 2016
 2017
 2018
 2019
 2020

Convert to DataFrame

Converting the vector of Pseudodata to DataFrame:

df = convert(DataFrame, pdata)
first(df,5)

Compute the median and the 1-α confidence interval for each of the configurations

α = 0.5

df[:, :y] = median.(df.Distribution)
df[:, :ymin] = quantile.(df.Distribution, α/2)
df[:, :ymax] = quantile.(df.Distribution, 1-α/2)

first(df, 5)

Display the discharge of all configurations

set_default_plot_size(12cm, 8cm)
plot(df, x=:Year, y=:y, color=:Configuration, Geom.line,
    ymin=:ymin, ymax=:ymax, Geom.ribbon,
    Guide.ylabel("Discharge [m³/s]"))

Display the discharge of one specific configuration

config = "MG24HQ"

df2 = filter(row -> row.Configuration ==config, df)

set_default_plot_size(12cm, 8cm)
plot(df2, x=:Year, y=:y, Geom.line, Geom.point,
    ymin=:ymin, ymax=:ymax, Geom.ribbon,
    Guide.ylabel("Discharge [m³/s]"))

Stationary extreme value model

Statistical model definition

model1 = PseudoMaximaModel(pdata, prior=[Flat(), Flat(), LocationScale(-.5, 1, Beta(6,9))])

PseudoMaximaModel
  datadistribution: Vector{Pseudodata}[6]
  location: μ ~ 1
  logscale: ϕ ~ 1
  shape: ξ ~ 1
  prior: [Mamba.Flat(), Mamba.Flat(), Distributions.LocationScale{Float64, Distributions.Beta{Float64}}(
μ: -0.5
σ: 1.0
ρ: Distributions.Beta{Float64}(α=6.0, β=9.0)
)
]

Sampling of the posterior distribution by MCMC

fm1 = fitbayes(model1)

PseudoMaximaEVA
   model: PseudoMaximaModel
   maxima:
		Mamba.Chains
		Iterations :		10000:20000
		Thinning interval :	10
		Chains :		1
		Samples per chain :	1001
		Value :			Array{Float64, 3}[1001,60,1]
   parameters:
		Mamba.Chains
		Iterations :		10000:20000
		Thinning interval :	10
		Chains :		1
		Samples per chain :	1001
		Value :			Array{Float64, 3}[1001,3,1]

Display of the annual maximum estimates

# Compute the estimates
ŷ = vec(mean( fm1.maxima[:, :, 1].value, dims=1))

# Compute the 95% credible interval
ymin = [quantile(vec(fm1.maxima[:, j, 1].value), .025) for j in 1:length(ŷ)]
ymax= [quantile(vec(fm1.maxima[:, j, 1].value), .975) for j in 1:length(ŷ)]

df = DataFrame(Year = years, Discharge=ŷ, ymin = ymin, ymax=ymax)

# Plot the ponctual estimates and the intervals
set_default_plot_size(12cm, 8cm)
plot(df, x=:Year, y=:Discharge, Geom.line, Geom.point,
    ymin=ymin, ymax=ymax, Geom.ribbon
)

Fit of the GEV law adjusted to the estimates of the annual maxima

General fit with residuals

set_default_plot_size(21cm ,16cm)
diagnosticplots(fm1, step = 20)